      SUBROUTINE  CPRSWF(C,XI,M,N,MINT,INDX,NMAX,MAXJW,JW,D,ALAMD,      00010000
     +                   SIGMA,DSUM,R1,DR1,R2,DR2)                      00020000
C                                                                       00030000
CC        THE PROLATE SPHEROIDAL RADIAL FUNCTIONS WITH COMPLEX ARGUMENTS00040000
CC        THE SIZE PARAMETER C IS COMPLEX.                              00050000
CC        INDX=1,  EXTERNAL FIELD*C1   ;   INDX=2,  INTERNAL FIELD*C2   00060000
CC        MINT = INDEX PARAMETER FOR CALCULATION OF BESSEL FUNCTIONS.   00070000
CC        M .LE. 40     :     JW .LE. 150                               00080009
CC        THE SECOND KIND FUNCTIONS ARE NEVER CALCULATED FOR INDX = 2.  00090000
C       ------------------------------  REVISED IN 1978 --------------  00100000
      PARAMETER  (JMX=150,MXPQ=190,MXSJ=150,MXSN=190)                   00110007
      COMPLEX*16 DX,W,T1,T2,T3,RATIO,SIGMA,PI,ALAMD,DSUM,E,F,G,R1A,DR1A,00120000
     1           RCONV1,R1,DR1,R2A,DR2A,RCONV2,R2N,DR2N,DR2NC,QSUM,DQSUM00130000
     2          ,QD,DN0,DP0,COEF2,PSUM,DPSUM,R2Q,DR2Q,DR2QC,AINT,C,C2,  00140000
     3           R2B,DR2B,RCONV3,TN,DTN,TN2,DTN2,CF,DCF,SUM,SUMD,V,DV,  00150003
     4           R2,DR2,DN(80),DNRAT(80),DP(200),DPRAT(200),D(JMX),     00160007
     5           SJ(MXSJ,2),DSJ(MXSJ,2),SN(MXSN,2),DSN(MXSN,2),AMX(JMX) 00170005
      REAL*8  Q,QNA,QNB, DQN,ALEG,XI,XI2,XIA,XIB,SJFACT,FLOAT,          00180012
     1        FF,PX(MXPQ),DPX(MXPQ),QX(MXPQ),DQX(MXPQ),QN(80),          00190007
     2        AL,ALD,B,BD,S1,S2,S3,RFAC1,RFAC2,RFAC3,RFAC4,FACTMM       00200011
      FLOAT(N)=DFLOAT(N)                                                00210000
      E(J)=C2*DFLOAT((M+M+J)*(M+M+J-1))/DFLOAT((2*(M+J)-1)*(2*(M+J)+1)) 00220000
      F(J)=DFLOAT((M+J)*(M+J+1))+C2*DFLOAT(2*(M+J)*(M+J+1)-2*M*M-1)/    00230000
     +     DFLOAT((2*(M+J)-1)*(2*(M+J)+3))-ALAMD                        00240000
      G(J)=C2*DFLOAT((J+1)*(J+2))/DFLOAT((2*(M+J)+1)*(2*(M+J)+3))       00250000
C                                                                       00260000
      DR2NC=(0.D0,0.D0)                                                 00270000
      DR2QC=(0.D0,0.D0)                                                 00280000
      XI2=XI*XI-1.D0                                                    00290000
      C2=C*C                                                            00300000
      H=C                                                               00310000
      HXX=H*XI*XI                                                       00320000
      LMTN=0.835*HXX                                                    00330000
      L=N-M                                                             00340000
C                                                                       00350011
        IF(M.EQ.0)  THEN                                                00360000
      XIA=1.D0                                                          00370000
      XIB=0.D0                                                          00380000
        ELSE                                                            00390003
      XIA=DSQRT(XI2/(XI*XI))**M                                         00400000
      XIB=DFLOAT(M)/(XI*XI2)                                            00410000
        END IF                                                          00420000
       IF(MAXJW.GE.JMX-1) MAXJW= JMX-2                                  00430008
      R2=(0.D0,0.D0)                                                    00440000
      DR2=(0.D0,0.D0)                                                   00450000
      PSCRT=(1.0E-12)*(10.0**(L/10.0))                                  00460000
CC        SPHERUCAL BESSEL FUNCTIONS                                    00470000
      IF(M.EQ.MINT.AND.N.EQ.M) GO TO 4002                               00480000
      GO TO 4005                                                        00490000
 4002 NSJ=MAXJW+M                                                       00500000
       IF(NSJ.GE.MXSJ-1)  NSJ= MXSJ-2                                   00510010
      DX=C*XI                                                           00520000
      NUP=10+DX/5.                                                      00530000
      NSWT=NSJ                                                          00540000
      NA=NSJ+NUP                                                        00550000
      W=(CDSIN(DX)-DX*CDCOS(DX))/(DX*DX)                                00560000
      T3=(0.0D0,0.0D0)                                                  00570000
      T2=(1.D-70,1.D-70)                                                00580000
      IF(DABS(DIMAG(DX)).LT.1.0D-10) T2=(1.D-70,0.D0)                   00590000
      DO 4000 NB=1,NA                                                   00600000
      NC=NA+1-NB                                                        00610000
      T1=DFLOAT(NC+NC+1)*T2/DX-T3                                       00620000
      IF(CDABS(T1)-1.D70) 4001,4003,4003                                00630000
 4003 T1=T1*1.0D-70                                                     00640000
      T2=T2*1.0D-70                                                     00650000
      NSWT=NC                                                           00660000
 4001 T3=T2                                                             00670000
      T2=T1                                                             00680000
      IF(NC.GT.NSJ) GO TO 4000                                          00690000
      SJ(NC,INDX)=T1                                                    00700000
 4000 CONTINUE                                                          00710000
      RATIO=W/SJ(2,INDX)                                                00720000
      DO 4010 ND=1,NSJ                                                  00730000
      IF(ND.GT.NSWT) GO TO 4015                                         00740000
      SJ(ND,INDX)=RATIO*SJ(ND,INDX)                                     00750000
      GO TO 4010                                                        00760000
 4015 SJFACT=1.D-70                                                     00770000
      IF(CDABS(SJ(ND,INDX)).LT.1.0D-5) SJFACT=0.D0                      00780000
      SJ(ND,INDX)=RATIO*SJ(ND,INDX)*SJFACT                              00790000
 4010 CONTINUE                                                          00800000
      DSJ(1,INDX)=-SJ(2,INDX)                                           00810000
      DO 4020 NE=2,NSJ                                                  00820000
      NF=NE-1                                                           00830000
 4020 DSJ(NF+1,INDX)=SJ(NF,INDX)-DFLOAT(NF+1)*SJ(NF+1,INDX)/DX          00840000
C         SPHERICAL NEUMANN FUNCTIONS                                   00850000
      IF(INDX.EQ.2) GO TO 4005                                          00860000
      IF(XI.LT.1.005) GO TO 4005                                        00870000
      IF(HXX.LT.10.0) GO TO 4005                                        00880000
      SCLSN=1.0E-20                                                     00890000
      IF(CDABS(DX).GE.15.0) SCLSN=1.0E-10                               00900000
      SN(1,INDX)=-CDCOS(DX)/DX                                          00910000
      SN(2,INDX)=-(CDCOS(DX)+DX*CDSIN(DX))/(DX*DX)                      00920000
      SN(1,INDX)=SN(1,INDX)*SCLSN                                       00930000
      SN(2,INDX)=SN(2,INDX)*SCLSN                                       00940000
      DSN(1,INDX)=-SN(2,INDX)                                           00950000
      MAXSN=MAXJW+M                                                     00960000
       IF(MAXSN.GE.MXSN-1) MAXSN= MXSN-2                                00970007
      DO 4030 NP=2,MAXSN                                                00980000
      NQ=NP-1                                                           00990000
      IF(NQ.LE.1) GO TO 4040                                            01000000
      SN(NQ+1,INDX)=DFLOAT(NQ+NQ-1)*SN(NQ,INDX)/DX-SN(NQ-1,INDX)        01010000
 4040 DSN(NQ+1,INDX)=SN(NQ,INDX)-DFLOAT(NQ+1)*SN(NQ+1,INDX)/DX          01020000
 4030 CONTINUE                                                          01030000
 4005 CONTINUE                                                          01040000
CC        ASSOCIATED LEGNDRE FUNCTIONS                                  01050000
      IF(INDX.EQ.2) GO TO 4007                                          01060000
      IF(N.GT.M) GO TO 4007                                             01070000
      MAXPQ=MAXJW+M                                                     01080000
       IF(MAXPQ.GE.MXPQ-1) MAXPQ= MXPQ-2                                01090007
      MLEG=0                                                            01100000
      ALEG=1.D0                                                         01110000
 4070 MLEG=MLEG+1                                                       01120000
      ALEG=ALEG*DFLOAT(2*MLEG-1)                                        01130000
      IF(MLEG.LT.M)  GO TO 4070                                         01140000
      PX(1)=ALEG*(DSQRT(XI2)**M)                                        01150000
      PX(2)=DFLOAT(2*M+1)*XI*PX(1)                                      01160000
          CALL  ALEGQ(XI,M,MAXPQ,QX)                                    01170000
      DO 4050 I=1,MAXPQ                                                 01180000
      MI=M+I                                                            01190000
      PX(I+2)=(DFLOAT(MI+MI+1)*XI*PX(I+1)-DFLOAT(M+MI)*PX(I))/FLOAT(I+1)01200000
      DPX(I)=(DFLOAT(I)*PX(I+1)-DFLOAT(MI)*XI*PX(I))/XI2                01210000
      IF(I.EQ.1) GO TO 4050                                             01220000
      DQX(I-1)=(DFLOAT(I-1-M)*QX(I)-DFLOAT(I-1)*XI*QX(I-1))/XI2         01230000
 4050 CONTINUE                                                          01240000
      IF(M.EQ.0) GO TO 4007                                             01250000
      S1=QX(2)                                                          01260000
      S2=QX(1)                                                          01270000
      DO 4080 KA=1,M                                                    01280000
      KK=KA-1                                                           01290000
      S3=(DFLOAT(1-2*KK)*XI*S2+DFLOAT(KK+M-1)*S1)/DFLOAT(M-KK)          01300000
      S1=S2                                                             01310000
      S2=S3                                                             01320000
      QN(KA)=S3                                                         01330000
 4080 CONTINUE                                                          01340000
 4007 CONTINUE                                                          01350000
CC        EIGENVALUES AND EXPANSION COEFFICIENTS (D) - HODGE METHOD     01360000
      IF(N.GT.M) GO TO 4009                                             01370000
         CALL  CHODGE(M,C,NMAX,AMX, 1)                                  01380000
 4009 MTRX=L+1                                                          01390000
      AINT=AMX(MTRX)                                                    01400000
C                                                                       01410010
         CALL  CDCOEF(C,M,N,1,AINT,JW,D,ALAMD,PI,SIGMA,DSUM,FCTR)       01420010
      IF(CDABS(DSUM).LT.1.0D-50) GO TO 3                                01430000
CC            RADIAL FUNCTION AND ITS DERIVATIVE - RJ(M,N,HX),DRJ       01440000
CC            THE RADIAL FUNCTIONS OF FIRST KIND                        01450000
      IFIN=JW                                                           01460000
      IF(MOD(L,2).EQ.0) IFIN=JW+1                                       01470000
      R1A=(0.D0,0.D0)                                                   01480000
      DR1A=(0.D0,0.D0)                                                  01490000
      DO 1060 JR=1,IFIN,2                                               01500000
      J=JR-1                                                            01510000
      IF(MOD(L,2).EQ.1) J=JR                                            01520000
      MR=M+J                                                            01530000
      IF(MR.GT.NSJ) GO TO 1060                                          01540000
      IR=IABS(J-L)                                                      01550000
      Q=-1.D0                                                           01560000
      IF(MOD(IR,4).EQ.0) Q=1.D0                                         01570000
      RFAC1= 1.D0                                                       01580013
      RFAC2= 1.D0                                                       01590010
        IF(M.GE.1)  THEN                                                01600010
      DO 15 IP=1, M                                                     01610010
      RFAC1=RFAC1*DFLOAT(J+IP)                                          01620010
   15 RFAC2=RFAC2*DFLOAT(J+IP+M)                                        01630010
        END IF                                                          01640010
      R1A=R1A+Q*D(JR)*RFAC1*SJ(MR+1,INDX)*RFAC2/FCTR                    01660013
      RCONV1= Q*D(JR)*RFAC1*DSJ(MR+1,INDX)*RFAC2/FCTR                   01661013
      DR1A=DR1A+RCONV1                                                  01670011
       IF(CDABS(RCONV1/DR1A).LE.1.0D-15) GO TO 5                        01680010
 1060 CONTINUE                                                          01690000
    5 R1= XIA*R1A/PI                                                    01700013
      DR1= XIB*R1+C*XIA*DR1A/PI                                         01710013
      IF(INDX.EQ.2) GO TO 4                                             01720000
C           RADIAL FUNCTION OF SECOND KIND IN TERMS OF SPHERICAL NEUMAN 01730000
C           SERIES WITH INTEGRAL METHOD BY SINHA & MACPHIE(1974)        01740000
      DFN=1.D+10                                                        01750000
      R2N= 0.D0                                                         01760000
      DR2N= 0.D0                                                        01770000
      IF(XI.LT.1.005) GO TO 7                                           01780000
      IF(HXX.LT.10.0) GO TO 7                                           01790000
      IF(L.GT.LMTN) GO TO 7                                             01800000
      IFN=IFIN                                                          01810000
    6 IFN=IFN-2                                                         01820000
      NR=IFN-2                                                          01830000
      R2A=0.D0                                                          01840000
      DR2A= 0.D0                                                        01850000
      DO 1070 JR=1,IFN,2                                                01860000
      J=JR-1                                                            01870000
      IF(MOD(L,2).EQ.1) J=JR                                            01880000
      MR=M+J                                                            01890000
      IF(MR.GE.MAXSN) GO TO 6                                           01900000
      IR=IABS(J-L)                                                      01910000
      Q=-1.D0                                                           01920000
      IF(MOD(IR,4).EQ.0) Q=1.D0                                         01930000
      RFAC3= 1.D0                                                       01940000
      RFAC4= 1.D0                                                       01950010
        IF(M.GE.1)  THEN                                                01960010
      DO 17 IP=1, M                                                     01970010
      RFAC3=RFAC3*DFLOAT(J+IP)                                          01980010
   17 RFAC4=RFAC4*DFLOAT(J+IP+M)                                        01990010
        END IF                                                          02000010
      RCONV2= Q*D(JR)*RFAC3*SN(MR+1,INDX)*RFAC4                         02010014
      RCONV3= Q*D(JR)*RFAC3*DSN(MR+1,INDX)*RFAC4                        02020014
      IF (JR.EQ.NR) TN=  RCONV2                                         02030000
      IF (JR.EQ.NR) DTN=RCONV3                                          02040000
      IF (JR.EQ.IFN) TN2=  RCONV2                                       02050000
      IF (JR.EQ.IFN) DTN2=RCONV3                                        02060000
      IF (JR.EQ.IFN) GO TO 8                                            02070000
      R2A=R2A+RCONV2                                                    02080000
      DR2A=DR2A+RCONV3                                                  02090000
 1070 CONTINUE                                                          02100000
    8 IF (MOD(L,2).EQ.0) NR=NR-1                                        02110000
         CALL CEM(M,NR,C,XI,ALAMD,V,AL,SUM,B,0)                         02120000
         CALL CEM(M,NR,C,XI,ALAMD,DV,ALD,SUMD,BD,1)                     02130000
      CF=TN/DEXP(-AL*(M+NR))/DFLOAT(M+NR)**B/CDEXP(SUM)                 02140000
      DCF=DTN/DEXP(-ALD*(M+NR))/DFLOAT(M+NR)**BD/CDEXP(SUMD)            02150000
      R2B=(CF*DEXP(-AL*(M+NR+2))/AL*V+TN2)/2.                           02160000
      DR2B=(DCF*DEXP(-ALD*(M+NR+2))/ALD*DV+DTN2)/2.                     02170000
       R2N= XIA*(R2A+R2B)/PI/FCTR                                       02180010
       DR2N= XIB*R2N+C*XIA*(DR2A+DR2B)/PI                               02190010
      R2N=R2N/SCLSN                                                     02200000
      DR2N=DR2N/SCLSN                                                   02210000
      DR2NC=(1.0/(C*XI2)+R2N*DR1)/R1                                    02220000
      DFN=CDABS((DR2NC-DR2N)/DR2NC)                                     02230000
    7 CONTINUE                                                          02240000
CCC       SPHEROIDAL RADIAL FUNCTIONS OF THE SECOND KIND                02250000
CCC       IN TERMS OF THE ASSOCIATED LEGENDRE FUNCTIONS                 02260000
      IF(MOD(L,2).EQ.1) GO TO 3000                                      02270000
      JFIN=2*M                                                          02280000
      JINT=2                                                            02290000
      GO TO 3010                                                        02300000
 3000 JFIN=2*M-1                                                        02310000
      JINT=1                                                            02320000
 3010 QSUM=(0.D0,0.D0)                                                  02330000
      DQSUM=(0.D0,0.D0)                                                 02340000
      DO 10 I=1,IFIN,2                                                  02350000
      IA=I-1                                                            02360000
      IF(MOD(L,2).EQ.1)  IA=I                                           02370000
      MA=M+IA                                                           02380000
      IF(MA.GE.MAXPQ) GO TO 11                                          02390000
      QD=D(I)*DQX(MA+1)                                                 02400000
      QSUM=QSUM+D(I)*QX(MA+1)                                           02410000
      DQSUM=DQSUM+QD                                                    02420000
      IF(CDABS(QD/DQSUM).LT.1.0D-13) GO TO 11                           02430000
   10 CONTINUE                                                          02440000
   11 IF(M.EQ.0) GO TO 3030                                             02450000
      DNRAT(JFIN+2)=(0.D0,0.D0)                                         02460000
      DO 20 J=JINT,JFIN,2                                               02470000
      JA=JINT+JFIN-J                                                    02480000
   20 DNRAT(JA)=-E(-JA+2)/(F(-JA)+G(-JA-2)*DNRAT(JA+2))                 02490000
      DN0=D(1)                                                          02500000
      DO 30 JB=JINT,JFIN,2                                              02510000
      DN(JB)=DN0*DNRAT(JB)                                              02520000
      MJB=M-JB                                                          02530000
      MJC=MJB+1                                                         02540000
      IF(MJB.LT.0) GO TO 2050                                           02550000
      QNA=QX(MJB+1)                                                     02560000
      QNB=QX(MJC+1)                                                     02570000
      GO TO 2060                                                        02580000
 2050 MJD=-MJB                                                          02590000
      IF(MJD.EQ.1) GO TO 2070                                           02600000
      QNA=QN(MJD)                                                       02610000
      QNB=QN(MJD-1)                                                     02620000
      GO TO 2060                                                        02630000
 2070 QNA=QN(1)                                                         02640000
      QNB=QX(1)                                                         02650000
 2060 DQN=(DFLOAT(MJB-M+1)*QNB-DFLOAT(MJB+1)*XI*QNA)/XI2                02660000
      QSUM=QSUM+DN(JB)*QNA                                              02670000
      DQSUM=DQSUM+DN(JB)*DQN                                            02680000
      DN0=DN(JB)                                                        02690000
   30 CONTINUE                                                          02700000
      IF(MOD(L,2).EQ.1) GO TO 3040                                      02710000
      COEF2=FF(M,N)*FACTMM(M)*DN(JFIN)*C*PI/(DFLOAT(M+M-1)*C**M)*       02720010
     +      FCTR                                                        02730010
      GO TO 3050                                                        02740000
 3040 COEF2=-FF(M,N)*FACTMM(M)*DN(JFIN)*C2*PI/(DFLOAT((M+M-3)*(M+M-1))* 02750010
     +      C**M)*FCTR                                                  02760010
 3050 CONTINUE                                                          02770000
      GO TO 3060                                                        02780000
 3030 IF(MOD(L,2).EQ.1) GO TO 3070                                      02790000
      COEF2=-FF(M,N)*D(1)*C*PI                                          02800000
      GO TO 3060                                                        02810000
 3070 COEF2=-FF(M,N)*D(1)*C2*PI/3.0D0                                   02820000
 3060 CONTINUE                                                          02830000
      IH=H                                                              02840000
      IF(MOD(L,2).EQ.1) GO TO 3080                                      02850000
      KINT=M+M+2                                                        02860000
      KFIN=30+2*(M+IH)                                                  02870000
      GO TO 3090                                                        02880000
 3080 KINT=M+M+1                                                        02890000
      KFIN=31+2*(M+IH)                                                  02900000
 3090 IF(KFIN.GE.138) KFIN=138+MOD(L,2)                                 02910007
 2040 DPRAT(KFIN+2)=(0.D0,0.D0)                                         02920000
      DP0=D(1)                                                          02930000
      IF(M.NE.0)  DP0=DN(JFIN)                                          02940000
      DO 40 KA=KINT,KFIN,2                                              02950000
      KB=KINT+KFIN-KA                                                   02960000
      IF(KB.EQ.KINT) GO TO 2010                                         02970000
      DPRAT(KB)=-E(-KB+2)/(F(-KB)+G(-KB-2)*DPRAT(KB+2))                 02980000
      GO TO 40                                                          02990000
 2010 IF(MOD(L,2).EQ.1) GO TO 2020                                      03000000
      DPRAT(KB)=C2/FLOAT((M+M+1)*(M+M-1))/(F(-KB)+G(-KB-2)*DPRAT(KB+2)) 03010000
      GO TO 40                                                          03020000
 2020 DPRAT(KB)=-C2/FLOAT((M+M-1)*(M+M-3))/(F(-KB)+G(-KB-2)*DPRAT(KB+2))03030000
   40 CONTINUE                                                          03040000
      PSUM=(0.D0,0.D0)                                                  03050000
      DPSUM=(0.D0,0.D0)                                                 03060000
      DO 50 KC=KINT,KFIN,2                                              03070000
      DP(KC)=DP0*DPRAT(KC)                                              03080000
      KD=KC-M-1                                                         03090000
      KE=KD+1-M                                                         03100000
      IF(KE.GT.MAXPQ) GO TO 50                                          03110000
      PSUM=PSUM+DP(KC)*PX(KE)                                           03120000
      DPSUM=DPSUM+DP(KC)*DPX(KE)                                        03130000
      DP0=DP(KC)                                                        03140000
   50 CONTINUE                                                          03150000
      KFINM=KFIN-M-M                                                    03160000
      IF(KFINM.GE.MAXPQ) GO TO 2030                                     03170000
      IF(KFIN.GE.137) GO TO 2030                                        03180007
      CRIT=CDABS(DP(KFIN)*PX(KFINM)/PSUM)                               03190000
      KFIN=KFIN+4                                                       03200000
      IF(CRIT.GT.PSCRT) GO TO 2040                                      03210000
 2030 KFIN=KFIN-4                                                       03220000
      R2Q=(QSUM+PSUM)/COEF2                                             03230000
      DR2Q=(DQSUM+DPSUM)/COEF2                                          03240000
      DR2QC=(1.0/(C*XI2)+R2Q*DR1)/R1                                    03250000
      DFQ=CDABS((DR2QC-DR2Q)/DR2QC)                                     03260000
      IF(XI.LT.1.005) GO TO 1                                           03270000
      IF(HXX.LT.10.0) GO TO 1                                           03280000
      IF(L.GT.LMTN) GO TO 1                                             03290000
      IF(DFQ.LE.DFN) GO TO 1                                            03300000
      R2=R2N                                                            03310000
      DR2=DR2N                                                          03320000
      GO TO 4                                                           03330000
    1 R2=R2Q                                                            03340000
      DR2=DR2Q                                                          03350000
      GO TO 4                                                           03360000
    3 WRITE(6,100)                                                      03370000
  100 FORMAT(1H0,20X,'*  CPRSWF  -  RETURN WITHOUT CALCULATION  *  COEFF03380000
     +ICIENTS D ARE ALL ZERO' )                                         03390000
    4 RETURN                                                            03400000
      END                                                               03410000
